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1. Introduction 

Intensity mapping of neutral hydrogen (HI) is a very promising, effieient tool to measure the large- 
seale strueture of the Universe (Battye et al. 2004; Vujanovie et al. 2009; Lidz et al. 2011; Peterson 
& Suarez 2012). The HI gas with speetral emission of A = 21em eaused by the spin-flip of the 
valenee eleetron is assumed to linearly traee the Dark Matter distribution. In intensity mapping, 
the entire HI flux is measured in large resolution elements whieh still allow observations of strue¬ 
ture on the B ary on Aeoustie Oseillation seale (BAO) (Wyithe et al. 2008; Chang et al. 2008; Bull 
et al. 2014a). These experiments are optimally eondueted by telescopes with a large field-of-view 
(FoV) through, for instance, multi-pixel feeds. The Green Bank telescope (GBT) has pioneered 
the detection of the large scale structure with the 21 cm signal (Chang et al. 2010; Masui et al. 
2013; Switzer et al. 2013). Near future experiments include single dish telescope such as BINGO 
(Battye et al. 2012) and interferometric designs such as CHIME (Bandura et al. 2014) or Tianlai 
(Chen 2012). The SKA will be able to carry out various intensity mapping experiments at different 
redshifts either in auto-correlation mode, i.e. single dish observations, or interferometric mode. 
For a more detailed discussion of intensity mapping with the SKA, we refer to the chapters Santos 
et al. (2014) and Bull et al. (2014b). 

The drawback of intensity mapping observations are the foregrounds mainly due to the high 
radio emission of our own Galaxy caused by synchrotron and free-free electron radiation. In the 
considered frequency range, the Galactic foregrounds significantly dominate the acquired maps by 
several orders of magnitude. The only full-sky observation of the Galactic radio emission is at 
frequency V = 408MHz (Haslam et al. 1982). This map combined with smaller surveys is the basis 
of most studies to understand and fully model the Galactic radio emission, i.e. de Oliveira-Costa 
et al. (2008); Jelic et al. (2008); Alonso et al. (2014). The varying spectral index with latitude 
and the lack of small-scale observations impede the comprehension of the Galactic foregrounds. 
The missing observational constraints on the Galactic foregrounds make the subtraction a crucial, 
highly sensitive step in the intensity mapping analysis pipeline. Foreground residuals can cause 
serious systematic effects which bias the cosmological analysis and therefore require careful con¬ 
sideration (Wolz et al. 2014). Most statistical methods employed for this task are based on the 
spectral smoothness of the Galactic foregrounds whereas the cosmological signal is expected to 
have low correlations between frequency bins. Extragalactic foregrounds caused by bright point 
source emission are expected to show different spatial structure due to gravitational clustering and 
be far less dominant than Galactic sources. Extragalactic contaminations are yet to be fully in¬ 
corporated in realistic simulations, however, not considered to cause any significant systematic 
effects. 

In this article, we review three different foreground subtraction approaches: the Karhunen- 
Eoeve transform (KE, see Shaw et al. 2014a,b), the singular value decomposition (SVD, see Masui 
et al. 2013) and the independent component analysis (FASTICA, see Wolz et al. 2014). The SVD 
and FASTICA are blind search methods which decompose the data in principal or independent 
components, respectively, whereas the KE transform separates foregrounds and 21cm signal by 
modeling their statistical properties. For the KE transform, we reconsider the foreground residuals 
of the power spectrum of an earlier study. We give an example of foreground removal by applying 
the FASTICA to an SKA phase 1 simulation for auto-correlation mode where we chose the settings 
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for SKAl-MID band 1. We emphasize that an SKAl-SUR band 2 experiment would have very 
similar prospects in terms of foreground removal and results are transferable to such an experiment. 
For SKA2 the foreground subtraction for experiments in the low redshift regime z ss 1 will be less 
challenging due to higher signal-to-noise ratio. Intensity mapping of structure in higher redshifts 
z PS 3, (which is of particular interest), is expected to be similarly challenging due to the increased 
sensitivity of the SKA2 experiment. We present the expected level of foreground residuals and 
BAO scale measurements in the angular power spectrum for FASTICA. In this work, we focus on 
methods which have been recently applied to intensity mapping simulations or data. In Chapman 
et al. (2014), foreground subtraction methods are presented in the context of SKA experiments on 
the epoch of reionization which potentially can be transferred to the intensity mapping framework. 

The article is structured as follows. In Sec. 2, we describe the 21 cm signal. Galactic fore¬ 
grounds and telescope noise properties including the simulation details. We also review the effect 
of instrumental errors on the foreground subtraction. In Sec. 3, the removal techniques are estab¬ 
lished and in the following Sec. 4, the performance of the FASTICA in foreground subtraction is 
evaluated by the considering the systematic errors. The conclusions are presented in Sec. 5. 

2. Description of the Signal Components 

In this section, we describe the nature of the cosmological signal we expect to observe and how we 
simulate it. We briefly outline the expected instrumental noise levels and review possible instru¬ 
mental errors. 

2.1 21cm Signal 

The intensity in a frequency bin 8v coming from the 21cm emission of an object at redshift z with 
neutral hydrogen mass Mhi, subtending a solid angle 8Q. is given by Abdalla & Rawlings (2005) 

^ 3hpAi2 1 Mhi .. 

\6nmn{{lFz)r{z)Y8v8Q^^^' 

where A 12 is the Einstein coefficient corresponding to the emission from the 21 cm hyperfine tran¬ 
sition, hp is Planck’s constant and is the hydrogen atom mass. Here, r(z) is the comoving 
curvature distance r(z) = csinn(//o x{z)/c)/{Hq Y^|f2j.|) and x{z) is the radial comoving 

distance x{z) = Io8z'/H{z'). 

This intensity /(v, n) can be written in terms of a black-body temperature in the Rayleigh-Jeans 
approximation T = Ic^/{IksV^), where kg is Boltzmann’s constant. Using this we can estimate 
the mean brightness temperature coming from redshift z and its fluctuations in terms of the neutral 
hydrogen density: 

T 21 (z, A) = (0.19055K) (1 + 5 hi) . (2.2) 

v/nM(l+z)3+nA 

Here vhi(z) is the neutral hydrogen fraction (i.e. fraction of the total baryon density in HI) and 5 hi 
is the HI overdensity field in redshift space (smoothed over the volume defined by 5v and 5fl). 

We simulated the cosmological signal by generating a three-dimensional realization of the neu¬ 
tral hydrogen density in the lightcone using a lognormal field as described in section 3 of Alonso 
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etal. (2014). This was done using generic ACDM cosmological parameters {Q.M,^b,^k,h,wo,Wa,o^,ns) 
(0.3,0.049,0,0.67,-1,0,0.8,0.96), and we further assumed a redshift dependence for the neutral 
hydrogen fraction of xhi = 0.008 (1 +z) and a clustering bias of 1 (5 hi = 5). Redshift-space dis¬ 
tortions were included using the radial velocity field inferred from the Gaussian overdensity field 
used for fhe lognormal realizafion. The simulafion box used is large enough fo encompass a full- 
sky volume fo redshiff 2.5 (wifh fhe observer placed af fhe cenfre), and has a spafial resolufion of 
abouf 2.65Mpc//j. 

We inferpolafed fhe density field info spherical femperature maps af differenl frequencies using 
Eq. 2.2. These maps were generated using fhe HEALPix package (Gorski el al. 2005) wifh a 
resolution parameter n_side = 512 (50 ~ 0.1 ty), al frequency intervals of 5v = 0.7 MHz befween 
405MHz (z ~ 2.5) and 945MHz (z ~ 0.5), making up a fofal of 770 frequency bands. 

2.2 Galactic Foregrounds 

Eor fhe simulafion used in Ihis analysis we generafed foreground realizalions for Galactic syn¬ 
chrotron emission and free-free emission (both Galactic and extragalactic). We used the method 
described in section 4 of Alonso et al. (2014) and did not include any leakage of the polarized 
synchrotron radiation. 

The free-free foregrounds were generated as a Gaussian random realization of the correspond¬ 
ing power spectra modelled in Santos et al. (2005). This models these foregrounds as isotropic 
processes, which is obviously not a good approximation for the Galactic case. However, due to 
their exceptionally smooth frequency dependence and subdominant amplitude, we do not believe a 
more sophisticated modelling is required at this stage. 

Galactic synchrotron is by far the largest foreground for intensity mapping, and a more com¬ 
plex method was used to simulate it. The method is largely based on that of Shaw et al. (2014b), 
and is also similar to those used in other studies (Jelic et al. 2010; Shaw et al. 2014a). The method 
starts by extrapolating the 408 MHz Haslam map (Haslam et al. 1982) to other frequencies using a 
given model for the direction-dependent synchrotron spectral index (for which we used the Planck 
Sky Model, see Delabrouille et al. 2013). Smaller-scale structure and frequency-decorrelation is 
then added on top of this through a Gaussian realization of the power spectra in Santos et al. (2005). 
Eurther details can be found in section 4.2 of Alonso et al. (2014). 


2.3 Instrumental Noise 


The noise RMS (flux sensifivify) for a single-poinfing observafion wifh a single-dish radio felescope 
is given by 


<ys = 


(2.3) 


IksTfiys 
Ag Svtp 

where Ag is fhe effective collecfing area of fhe dish, tp is fhe durafion of fhe observation, and we 
have assumed fhaf fhe noise is Gaussian and uncorrelafed. The fofal system femperafure is Tgys = 
Tinst + Tsky, where Tinst ~ few x lOK depends on fhe noise characteristics of fhe receiver sysfem, 
and Tsky ~ 60K(v/300MHz)^^-^^ accounfs for fhe femperafure of fhe sky due fo background radio 
emission. Converting Eq. (2.3) info a brighfness femperafure sensifivify gives Oj = Tsys/ y/5vtp in 
fhe Rayleigh-Jeans limit. 
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In our simulations, we assume that the various SKA eonfigurations will be able to perform 
10,000 hour auto-eorrelation surveys over an area of 30,000 sq. deg^ For an instrument with Adish 
independent dishes, the total survey area fltot can be eovered with an observation time per pointing 
of tp = ttot(Adish^B/^tot )5 where Q.b ~ is the beam solid angle of a single dish. Assuming 

no overlap between survey pointings, the brightness temperature sensitivity per pointing beeomes 


Oj 


T 

^ sys 

\/5vttot 


a 


'tot 


Adish^B 


(2.4) 


The noise RMS per pixel in our simulations is then given by apx = Oj / Opix 8 Vpix. 

We simulate an intensity mapping experiment whieh ean be performed by an SKAl-MID 
survey with band 1. We assume Tinst = 28K for Adish = 190 with a dish diameter of Ddish = 15m. 


2.4 Instrumental Effects 

In theory, the speetrally smooth foregrounds ean be very well extraeted by presented methods, how¬ 
ever, instrumental effeets significantly complicate the subtraetion. Most of the foreground subtrae- 
tion methods are not eapable of dealing with varying beams whieh is required for auto-eorrelation 
observations. The issue ean be mitigated by deeonvolving the data to the lowest resolution, as has 
been done for our simulation and for example in Wolz et al. (2014). Seale-dependent foreground 
removal methods are required for future experiments in order to not lose spatial information. 

The most prominent instrumental systematies are polarization leakage and frequeney-dependent 
beam distortions which cause spatial and speetral modes to mix. In Alonso et al. (2014), a first step 
is taken to inelude polarization leakage in the intensity mapping simulation such that future studies 
ean test the mode mixing effeets. Additionally, ealibration errors, telescope pointing errors and RFI 
adulterate signal proeessing. Existing methods need to be advaneed to inelude these instrumental 
effeets in the foreground removal step. 


3. Review of Foreground Removal Methods 

In the following seetion, we review seleeted foreground removal methods, the SVD, KL transform 
and FASTICA. Eaeh of these methods have been designed and applied to different types of intensity 
mapping data types and are part of power speetrum estimator pipelines. The KE transform is part 
of a detailed interferometer simulation for CHIME (Shaw et al. 2014a,b). The SVD is a more 
empirieal approaeh utilized to subtraet the foregrounds of GET observations (Masui et al. 2013; 
Switzer et al. 2013). The FASTICA has been tested on SKA-like simulations in Wolz et al. (2014). 
Eaeh of these methods work in a different mathematieal and experimental framework and we refer 
the reader to the respeetive artieles for more details on the baekground deseription. 

3.1 Karhunen-Loeve Transform 

The Karhunen-Eoeve transform has a long history in Cosmology (Bond 1995; Vogeley & Szalay 
1996), but was first suggested as a method for 21 em foreground removal in Shaw et al. (2014a). It 

*We focus on auto-correlation surveys here because of their better sensitivity to relevant cosmological scales for 
the SKA Phase 1 configurations; see Santos et al. (2014) for a discussion of the relative merits of auto-correlation and 
interferometric surveys. 
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has also been used for the analogous problem of E/B mode separation for polarisation of the CMB 
(Lewis et al. 2002; Bunn et al. 2003). In this section, we give an overview of the technique, and 
discuss its practical implementation in Shaw et al. (2014a,b). 

Ultimately foreground cleaning is simply a matter of finding a subset of our data within which 
there is significantly more 21 cm signal than there are astrophysical foregrounds. However, in the 
presence of mode mixing there is no immediately apparent representation which separates these 
signals. The Karhunen-Loeve transform (often called the Signal-to-Noise eigendecomposition), 
gives an automated way of deriving this basis from the two-point statistics of each component. 

This method requires models for the two-point statistics of both the signal and the foregrounds 
on the sky. We denote the matrix representation of the signal power spectrum as C 21 , whereas 
the foreground model, which includes both the synchrotron emission from our galaxy, and the 
contribution from a background of extragalactic point sources, is written as C/. Appropriate models 
are described in Shaw et al. (2014b) 

The Karhunen-Loeve transform seeks to find a linear fransformafion of fhe dafa d' = Pd such 
fhaf fhe covariance mafrices of fhe 21cm signal S = BC 2 iB^ and foregrounds F = BC/B^ are 
joinfly diagonalised. Thai is 

S^S' = PSP^=A, (3.1) 

and 

F ^ F' = PFPt = I, (3.2) 

where A is a diagonal mafrix, and I is fhe idenlily. In fhis diagonal basis we can simply compare 
fhe amounf of power expecfed in each mode by fhe rafio of fhe diagonal elemenfs (fhis is given 
by fhe corresponding enlries of A), and identify fhe regions of fhe space wilh low foreground 
confaminafion (large enlries in A). 

This fransformafion can be found by solving fhe generalised eigenvalue problem Sx = AFx. 
This gives a sef of eigenvectors x, and corresponding eigenvalues A. Wriling fhe eigenvecfors in a 
mafrix P, row-wise, gives fhe fransformafion mafrix fo diagonalise fhe covariances. The eigenval¬ 
ues X corresponding lo each eigenvector give fhe diagonal mafrix A. To isolale fhe 21 cm signal, we 
can simply selecf modes wilh eigenvalue (signal-fo-foreground power) greafer lhan some Ihresh- 
old. In Fig. 1 we show fhe effecl of foreground cleaning on fhe power speclrum errors of simulaled 
cylinder telescope similar lo Bandura el al. (2014). 

The Karhunen-Loeve Iransform is a general and very effective scheme for fhe componenl sep- 
aralion. Unforlunalely, fhe full covariance mafrices are dense, and very large, O(ApixAfreq) ~ 10^- 
10^ on a side, making fhe diagonalisalion, which is an 0{N^) operalion, impossible in fhe general 
case. However, fhis technique can still be applied to the restricted domain of transit telescopes, 
where the symmetries of the system allow the problem to be broken up into many smaller problems 
which are significantly more tractable. This decomposition is known as the m-mode formalism, and 
was first demonstrated in Shaw et al. (2014a). Overall the complexity is reduced to 
saving a factor of ~ 10^ in computation. 

3.2 Singular Value Decomposition 

In this section, we present the foreground cleaning formalism in Switzer et al. (2013) for the red- 
shifted HI survey with the Green Bank Telescope (GBT). It is based on the idea that the 21 cm sig- 


6 



Foreground Subtraction in Intensity Mapping 


Laura Wolz 



fcx / /iMpc ^ 


Figure 1: Forecast errors on the power spectrum as a fraction of its fiducial value for a 40m x 40m cylinder 
telescope similar to Bandura et al. (2014), observing between 400-500MHz. The three panels show the 
predicted errors without foregrounds (left), with unpolarised foregrounds (centre), and with fully polarised 
foregrounds (right). We can see that we lose sensitivity to the large scale line of sight mode. The dashed 
line indicates the predicted bound of the ‘foreground wedge’, showing that with perfect knowledge of our 
instrument, foregrounds can be successfully cleaned well into this region. 


nal is a line emission and has structure in the frequency/redshift domain, while the foregrounds are 
smooth on those scales. One can thus separate the two components in a model independent way by 
making use of principal component analysis. In this case, we arrange the three-dimensional maps 
into a frequency-frequency covariance matrix and use singular value decomposition to identify the 
dominant frequency modes. We label them as “foregrounds" and project out their contribution to 
obtain foreground cleaned maps. 

Specifically, we rearrange the three-dimensional map into an Ny x Nq matrix M, where Ng 
includes all two-dimensional spatial pixels. For the purpose of this comparison, we ignore thermal 
noise in the map. The empirical v — V covariance of the map is C = MM^/A^e, which contains 
both foregrounds and 21 cm signal. Making use of the component separation idea, the matrix can 
be factored as C = UAUT, where A is diagonal and contains the eigenvalues in descending order. 
We tag the first few modes as ’foregrounds’, and from each line of sight, we can then subtract a 
subset of the modes U that describe the largest components of the frequency variance through the 
operation (1 — USU^)M, where S is a selection matrix with 1 along the diagonal for modes to be 
removed and 0 elsewhere. 

In practice, the separation of foreground and signal modes is not perfect due to the imperfect 
characterization of instrumental response. The choice of the selection matrix S is a compromise 
between maximal foreground removal and minimal signal loss. To estimate the latter, we inject 
simulated 21cm signal to the data stream to determine the transfer function T which describes 
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loss of 21 cm signal due to foreground removal. As a rule of thumb, T = /’sig.out/Aig.in ~ [(1 — 
Am/A^v) (1 — Am/A-es)]^, where is the number of modes removed, Ny is the number of frequency 
channels and Ares is the number of angular resolution elements. A limited number of resolution 
elements can greatly reduce the efficacy of the foreground cleaning at the expense of signal. The 
details of foreground transfer function calculation and the signal compensation can be found in 
Switzer et al. (2013). 

Our approach to foreground removal is limited by the amount of information in the maps. The 
fundamental limitation here is not simply from the number of degrees of freedom along the line of 
sight, but instead is limited by the smaller of independent angular or frequency resolution elements 
in the map (Nityananda 2009). To see why this is the case, notice that in the absence of noise, 
our cleaning algorithm is equivalent to taking the SVD of the map directly: M = UEV^ and thus 
C oc MM^ = UE^U^, with the same set of frequency modes U appearing in both decompositions. 
The rank of C coincides with the rank of M and is limited by the number of either angular or 
frequency degrees of freedom. 


3.3 Independent Component Analysis 

In the following subsection, the basic principles of the fast independent component analysis FAS- 
TICA method (Hyvarinen 1999) are outlined. FASTICA has been successfully used for foreground 
removal in the astrophysical framework, see for example Maino et al. (2002); Bottino et al. (2010); 
Chapman et al. (2012). 

FASTICA is a method designed to decompose mixed signals into their independent constituents. 
Statistical independence is measured by the use of the Central Limit theorem which states that the 
probability distribution density (pdf) of a sum of independent variables is always more Gaussian 
than the pdf of one single component. The inverse application of that theorem implies that a 
single component can be extracted from the mixture by maximizing the non-Gaussianity of the 
pdf. This is turn implies that FASTICA can not be applied to Gaussian variables. We can formulate 
the problem in the linear equation x = As, where x is the signal, A the mixing matrix and s the 
independent components (IC). FASTICA blindly solves the inverted problem 


s = Wx (3.3) 

where the weighing matrix W defined as the inverse of A is unknown. As a measure of the non- 
Gaussianity of the ICs, FASTICA maximizes a proxy for the negentropy. 

FASTICA has been employed to subtract foregrounds in intensity mapping simulations in Wolz 
et al. (2014). Since the foregrounds are highly correlated between frequencies, FASTICA incorpo¬ 
rates them into the ICs. The 21 cm signal has a very low correlation which is close to Gaussian. 
Hence, FASTICA blindly reconstructs the foregrounds and the residuals of this analysis are cos¬ 
mological signal plus the receiver noise. In Wolz et al. (2014), it has been shown that foreground 
subtraction with FASTICA does not affect cosmological distance measures such as the BAG scale. 
An residual analysis of the SKA simulations with FASTICA is demonstrated in the following sec¬ 
tions. 
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4. Residual Errors and BAO recovery 

In this section, we present the foreground residuals and systematic errors when applying FASTICA 
to the described SKAl-MID simulation. 

4.1 Experimental Set-up 

The SKAl-MID band 1 simulation is based on the 21cm and Galactic foregrounds simulation 
presented in Alonso et al. (2014), where the frequency range is 406MHz < v < 945MHz with a 
frequency width of 5v = 0.7MHz per channel. The beam is approximated by a Gaussian beam with 
constant solid angle flbeam = 1.6deg^. The instrumental noise is simulated according to Equ. 2.4 
for an SKA-MID band 1 like observation. We mask out the sky above 63deg latitude in equatorial 
coordinates which results in a coverage of 30, OOOdeg^. The intensity maps are stacked into constant 
redshift bins with width 5z ~ 0.05 after the foreground removal to increase the signal-to-noise in 
the systematics analysis. 

4.2 Results 

We evaluate the foreground subtraction by estimating the angular power spectrum C{1) of the orig¬ 
inal 21 cm intensity maps and the reconstructed maps. We correct for the partial sky coverage with 
the Peebles approximation, as described in Wolz et al. (2014). In Fig. 2(a), we show the relative 
error CicA(^)/Corig(^) of the FASTiCA-cleaned maps for different number of ICs. The black error 
bars indicate the statistical variance of the measurements given by sampling error and instrumental 
noise. It can be seen that for a number of ICs higher than 4 the reconstruction converges towards 
the original input C{1) for scales i > 100. For smaller scales, the broadband spectrum is slightly 
distorted and we observe deviations of Cica(^) outside the statistical errors. 

In Fig. 2(b), the relative systematic error 5{£,Zi) is depicted for all redshift bins Zi on the 
x-axis and for multipole bins 5i = 20 in the y-direction. The systematic error is defined as 
-^sys = |Corig(^) — Cica(^)|- Wc divide that error by the statistical uncertainty on the measurements 
y^a{C{i)) composed of sampling variance and instrumental noise. The relative errors are pictured 
in a logarithmic scale such that values larger than zero imply insignificant systematic errors and 
smaller than zero indicate scales where systematics dominate the statistical errors. We see that the 
foregrounds at the edges of the redshift bins are poorly subtracted due to incomplete frequency 
information and those regions should be neglected in subsequent analysis. 

Furthermore, we present BAO scale measurements of the SKAl simulation in Fig. 3, the 
details of the methods can be found in Wolz et al. (2014). In Fig. 3(a), the BAO wiggles in the 
power spectrum produced by dividing by a BAO-featureless power spectrum model are shown for 
4 different redshift bins. The lines represent the fitted models to the estimated data points. The 
black measurements mark the simulation without any foreground contaminations and the red lines 
include systematic errors. We see that there is no significant shift of the BAO features induced 
by the foreground subtraction. Fig. 3(b) presents the scale distortion parameter a as a function 
of redshift where a = D^izi) /■ The BAO recovery is not significantly affected by the 
foreground subtraction for most redshifts Zi- 
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(a) Relative error on C(i) for 1.45 < z < 1.50 


(b) Relative systematic errors; FASTICA with IC=5 


Figure 2: The relative errors of the foreground subtraction with FASTICA on the angular power spectrum as 
a function of scale. On the left side, we present the fraction Cica{^)/C origi^) for different numbers of ICA 
for an example redshift bin. The relative systematic errors 5(£,z/) are given on a natural logarithmic scales 
in the right panel. 
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(a) BAO scale in angular power spectrum; FASTICA with IC=5 


(b) BAO fit; FASTICA with IC=5 


Figure 3: The BAO measurements of the original simulations without foregrounds (black) are compared to 
the measurements including foregrounds (red). In the left panel, we present the BAO wiggles in the angular 
power spectrum divided by a featureless fiducial model. The lines are produced by fitting a model to the 
measurements. The right panel shows the BAO scale distortion parameter a without and with foregrounds 
as a function of redshift. 
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5. Conclusions 

In this article, we have presented a review of three effective methods for subtracting foregrounds in 
intensity mapping experiments. We have simulated a future SKAl intensity mapping experiment 
with special focus on realistic foreground modeling. We presented a detailed residual analysis of 
the foreground-removed 21 cm signal with FASTICA. We conclude the article as follows. 

• Method comparison: The KL transform is based on a-priori knowledge on the proper¬ 
ties of foregrounds and 21 cm signal. SVD and FASTICA blindly decompose the observed 
data such that foregrounds are removed by subtracting an empirically chosen number of 
principal/independent modes. The blind methods are excellent in dealing with unexplored 
foregrounds or foreground fluctuations, however, prone to remove 21cm signal or leave 
foreground residuals in the data. By detailed, realistic simulations and residual analysis 
such as presented in this work, the errors can be significantly reduced and systematics well- 
understood. The KL transform has the advantage of separating the components based on 
their covariances which reduces described errors. This method is less versatile to exper¬ 
imental settings due to its high computational costs and can possibly be affected by poor 
foreground modeling. 

• Systematic errors: All presented methods have been shown to effectively clean intensity 
maps from high Galactic foregrounds in either simulations or observations such as the SVD 
method on the GBT data. Figs. 1 and 2 show that the systematic errors mostly affect the very 
large scales perpendicular to the line-of-sight. KL transform and FASTICA have been shown 
to not significantly bias the BAG detection in the cleaned data. We emphasize the importance 
of comprehensive and detailed simulations of intensity mapping experiments to ensure the 
reliability of the cosmological analysis. 

• Future requirements: Instrumental errors can influence the data analysis significantly by 
mode-mixing effects. The foreground removal can be coiTupted by telescope errors such as 
polarization leakage, telescope pointing and calibration errors. Shaw et al. (2014b) has in¬ 
cluded detailed instrumental effects in their simulation pipeline and studied the influence of 
each error on the power spectrum. Such effects require careful simulations for each individ¬ 
ual experimental set-up in order to understand the systematic effects and how they interact 
with the foreground removal. 

In this article, we outlined some of the latest efforts to subtract foregrounds from intensity mapping 
data. We showed that the 21 cm power spectrum can be successfully recovered via a wide range of 
scales and the BAO scale measurements are not biased by foreground subtraction. Within the next 
decade, those methods need to be advanced to the requirements of SKAl and SKA2 experiments 
with special attention to include polarization effects and instrumental errors into the systematic 
analysis. 
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